#### Figure 3 ####

#SSTOTAL graph
par(mar=c(5.1,4.1,1.1,3.5))
set.seed(10)
model<-sapply(seq(1,5,.1),function(x) rnorm(1,1+2*x))
plot(seq(1,5,.1),model,pch=16, xlab='x', ylab='y', type='n')
abline(h=mean(model), lty=2, col=1, lwd=2)
matguy<-cbind(seq(1,5,.1), mean(model), seq(1,5,.1), model)
apply(matguy,1,function(x) segments(x[1], x[2], x[3], x[4], lwd=2, col='dodgerblue'))
points(seq(1,5,.1), model, pch=16, cex=.9)
mtext(side=4,line=.5, expression(y==bar(y)), at=mean(model), cex=.9, col=1, las=2)
legend(1,12,expression(y[i]-bar(y)), lty=1, col='dodgerblue', cex=.95, lwd=2, bty='n')

#SSE graph
par(mar=c(5.1,4.1,1.1,5.1))
plot(seq(1,5,.1), model, pch=16, xlab='x', ylab='y', type='n')
lm(model~seq(1,5,.1))->out
abline(out,lty=2, col=1, lwd=2)
matguy<-cbind(seq(1,5,.1), fitted(out), seq(1,5,.1), model)
apply(matguy,1,function(x) segments(x[1], x[2], x[3], x[4], lwd=2, col='dodgerblue'))
points(seq(1,5,.1), model, pch=16, cex=.9)
mtext(side=4, line=.5, "regression\n line", at=coef(out)[1]+coef(out)[2]*5.2, cex=.8, col=1, las=2)
legend(1,12, expression(y[i]-hat(y)[i]), lty=1, col='dodgerblue', cex=.95, lwd=2, bty='n')

##### Figure 4 #####

#categorical predictors

H <- rep(0:1, c(20,20))
P <- rep(c(rep(0,10),rep(1,10)),2)
set.seed(20)
epsilon<-rnorm(40,0,1)
y<- 3+5*H+3*P+6*H*P+epsilon

lm(y~H+P+H:P)->out.h
summary(out.h)
lm(y~factor(H):factor(P)-1)->out.h1
summary(out.h1)
small.dat<-data.frame(H=c(1,2,1,2), P=c(1,1,2,2), mu=coef(out.h1))

plot(mu~P,data=small.dat,xlim=c(.8,2.2),axes=F,ylim=c(0,20), type='n', ylab=expression(paste('Mean response ',(mu))))
box()
axis(1,at=1:2,labels=c('0 (low)','1 (high)'))
axis(2)
cur.dat<-small.dat[small.dat$H==1,]
lines(cur.dat$P,cur.dat$mu,col=2,lwd=2)
cur.dat<-small.dat[small.dat$H==2,]
lines(cur.dat$P,cur.dat$mu,lwd=2,col=4)

arrows(.97,coef(out.h1)[1]+.3,.97,coef(out.h1)[2]-0.3,angle=45, length=.05, code=3,  col=1)
arrows(.97,0.3,.97,coef(out.h1)[1]-.3,angle=45, length=.05, code=3,  col=1)
segments(1,coef(out.h1)[1],2,coef(out.h1)[1],lty=2)
arrows(2.03,coef(out.h1)[3]+0.3,2.03,coef(out.h1)[2]+coef(out.h)[2]-.3,angle=45, length=.05, code=3,  col=1)
segments(1,coef(out.h1)[2],2,coef(out.h1)[2]+coef(out.h)[2],lty=2)
arrows(2.03,coef(out.h1)[1]+0.3,2.03,coef(out.h1)[3]-.3,angle=45, length=.05, code=3,  col=1)
arrows(2.03,coef(out.h1)[2]+coef(out.h)[2]+.3,2.03,coef(out.h1)[4]-0.3,angle=45, length=.05, code=3,  col=1)
arrows(2.03,0.3,2.03,coef(out.h1)[1]-0.3,angle=45, length=.05, code=3,  col=1)
text(.98,coef(out.h1)[1]/2,expression(beta[0]),pos=2,col=1)
text(.98,coef(out.h1)[1]+coef(out.h)[2]/2,expression(beta[1]),pos=2,col=1)
text(2.02,coef(out.h1)[1]/2,expression(beta[0]),pos=4,col=1)
text(2.02,coef(out.h1)[1]+coef(out.h)[3]/2,expression(beta[2]),pos=4,col=1)
text(2.02,coef(out.h1)[2]+coef(out.h)[3]/2,expression(beta[1]),pos=4,col=1)
text(2.02,coef(out.h1)[2]+coef(out.h)[3]+coef(out.h)[4]/2,expression(beta[3]),pos=4,col=1)

mycol=c(2,4)
mypch=c(15,16)
points(small.dat$P,small.dat$mu,pch=mypch[small.dat$H],col=mycol[small.dat$H])
legend('topleft',c('0 (low)','1 (high)'), col=mycol, pch=mypch, pt.cex=1, cex=.9,bty='n', title='H')
abline(h=0,lty=2)

mtext(expression(mu==beta[0]+beta[1]*H+beta[2]*P+beta[3]*group('(',H%*%P,')')),side=3,line=.2)